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Nonlinear  Corrections  to  Temperature  in  Computer  Simulations  of  Complex  Systems 


Abstract 

Temperature  is  the  familiar  thermodynamic  quantity  that  governs  heat  flow  between  large 
systems.  However,  temperature  comes  from  just  the  linear  (first-order)  derivative  of  entropy  with 
respect  to  energy,  so  that  nonlinear  corrections  may  contribute  significantly  to  the  equilibrium 
properties  of  small  systems.  Moreover,  the  nonlinear  corrections  also  influence  the  thermal  and 
dynamic  properties  of  independently-relaxing  nanometer-sized  regions  (“hot  spots”)  inside  bulk 
materials.  Several  experimental  techniques  have  shown  that  such  localized  regions  dominate  the 
primary  response  of  most  materials.  During  the  grant  period  we  have  greatly  extended  our 
fundamental  understanding  of  these  regions,  and  can  now  simulate  several  of  their  properties. 
This  Final  Report  describes  our  key  findings,  with  an  emphasis  on  comparing  computer 
simulations  to  experimental  data. 

Foreword 

Most  of  the  results  presented  here  come  from  Monte  Carlo  (MC)  simulations  of  the  Ising 
model.  The  first  nonlinear  correction  we  found  allows  simulations  that  exhibit  extensive  thermal 
fluctuations,  consistent  with  measured  heat  capacity  and  usual  definitions  of  extensive  entropy. 
Furthermore,  this  correction  significantly  improves  the  agreement  between  the  Ising  model  and 
measured  non-classical  scaling  behavior  of  ferromagnetic  materials  and  critical  fluids.  A  related 
nonlinear  correction  yields  spectra  that  exhibit  1/f-like  noise,  providing  a  fundamental 
mechanism  for  the  low-frequency  thermal  fluctuations  exhibited  by  most  materials.  Indeed  these 
corrections  give  a  unified  picture  for  several  empirical  laws,  including  commonly  measured 
deviations  from  the  mathematical  formulas.  A  simplified  Creutz  model  is  found  that  requires  two 
types  of  nonlinear  corrections  to  temperature.  The  model  consists  of  a  spin  system  coupled  to  an 
explicit  heat  bath.  One  mismatch  in  temperature  comes  from  finite-size  effects  in  the  heat  bath, 
while  the  second  mismatch  arises  if  excitations  in  the  heat  bath  are  localized  on  an  intermediate 
length.  Results  from  this  model  give  good  agreement  with  the  excess  specific  heat  of  imperfect 
crystals  at  low  temperatures.  Another  model,  which  involves  adiabatic  demagnetization  with 
weak  coupling  to  the  heat  bath,  describes  measurements  of  time-dependent  magnetoresistance  in 
nanowires  at  low  temperatures.  Recently  we  have  found  that  fluctuations  in  molecular-dynamics 
(MD)  simulations  of  dilute  systems  at  low  temperatures  deviate  from  the  equipartition  theorem. 
We  are  currently  searching  for  the  optimal  nonlinear  corrections  for  these  MD  simulations,  so 
that  they  may  also  yield  improved  agreement  with  thermodynamic  behavior  and  with 
experiments,  similar  to  the  corrections  we  have  found  for  MC  simulations.  The  fundamental  goal 
of  our  research  is  to  fully  understand  finite-size  thermal  effects  inside  complex  materials,  with  a 
practical  result  of  improved  accuracy,  efficiency,  and  predictive  power  of  computer  simulations. 
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1.  Introduction  [1] 

Temperature  is  the  primary  thermodynamic  quantity  that  governs  the  equilibrium  of  systems 
in  thermal  contact  with  a  heat  bath.  Ideally,  the  heat  bath  must  be  homogeneous  and  effectively 
infinite,  while  the  coupling  to  the  system  must  be  weak,  so  that  the  bath  changes  only  from  heat 
exchange  with  the  system.  We  have  investigated  some  consequences  of  small  thermal  baths,  with 
strong  coupling  to  the  system.  We  are  guided  by  the  theory  of  small-system  thermodynamics, 
introduced  by  Hill  in  1962  [2,3],  which  establishes  the  fundamental  laws  that  govern  finite-size 
effects  in  thermodynamics.  Although  originally  developed  to  describe  individual  molecules  and 
isolated  nanoparticles,  we  use  this  “nanothermodynamics”  for  treating  the  heterogeneous 
distribution  of  independently  relaxing  regions  inside  bulk  materials.  Several  experimental 
techniques  have  shown  that  such  nanometer-sized  regions  dominate  the  primary  response  of 
most  materials  [4-9]. 

We  use  two  key  concepts  from  nanothermodynamics:  Hill’s  subdivision  potential  (E)  and  the 
fully-open  nanocanonical  ensemble.  The  nanocanonical  ensemble  facilitates  systematic  treatment 
of  independently-relaxing  regions  in  thermal  contact  with  an  ensemble  of  similar  regions, 
consistent  with  the  primary  response  of  most  materials.  Other  ensembles  impose  artificial 
internal  constraints  on  the  regions,  as  shown  schematically  in  Fig.  (1).  In  general,  the  distribution 
of  internal  region  sizes  is  governed  by  Hill’s  E.  Furthermore,  this  E  must  be  included  in  the  1st 
law  of  thermodynamics  if  total  energy  is  to  be  conserved  [10].  Hill’s  E  can  be  understood  by 
comparison  to  Gibbs’  chemical  potential,  /jl.  //  is  the  change  in  energy  to  take  a  single  particle 
from  a  bath  of  particles  into  the  system,  whereas  E  is  the  change  in  energy  to  take  a  cluster  of  n 
interacting  particles  from  a  bath  of  clusters  into  the  system,  and  in  general  n  interacting  particles 
do  not  have  the  same  energy  as  n  isolated  particles  due  to  surface  terms,  length-scale  effects, 
thermal  fluctuations,  etc.  Thus,  E  is  needed  to  systematically  treat  the  nonlinear  and  non- 
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Figure  1.  From  Ref.  [12].  Various  ensembles  for  fluctuations 
inside  a  bulk  sample.  The  microcanonical  ensemble  applies 
to  fully-closed  fluctuations  that  conserve  number  of  particles 
(n),  volume  (v),  and  energy  (e).  The  canonical  ensemble 
(, n,v,T)  applies  to  fluctuations  at  constant  volume  when  only 
heat  flows  in  and  out  from  the  rest  of  the  sample.  The  grand 
canonical  ensemble  (ju,v,T)  applies  to  fluctuations  at  constant 
volume  when  particles  also  flow  in  and  out  from  the  rest  of 
the  sample.  The  “nanocanonical”  ensemble  (/ i,P,T)  applies  to 
fully-open  fluctuations,  where  the  volume  of  the  fluctuation 
is  allowed  to  change  as  particles  flow  in  and  out. 


4 


extensive  contributions  to  energy  from  a  system  that  contains  a  heterogeneous  distribution  of 
independent  regions.  Moreover,  £-> 0  governs  the  distribution  of  regions  if  their  sizes  are  not 
fixed,  similar  to  how  //->()  governs  the  distribution  of  particles  if  their  numbers  are  not  fixed. 

We  have  used  the  concepts  of  nanothermodynamics  as  a  guide  to  treat  independent  thermal 
fluctuations  inside  bulk  samples.  We  find  that  the  resulting  models  and  theories  yield  improved 
agreement  with  the  measured  response  of  many  materials  [11-17].  Indeed,  these  concepts 
provide  a  physical  explanation  for  several  empirical  formulas  that  are  used  to  characterize  the 
complex  dynamics  of  realistic  materials.  In  fact  we  have  shown  that  nanothermodynamics 
provides  a  fundamental  basis  for:  non-exponential  relaxation,  non-Arrhenius  activation,  non- 
classical  critical  scaling,  non-homogeneous  response,  and  non-Nyquist  noise,  including  most- 
often  observed  deviations  from  the  simple  empirical  formulas.  First  I  will  describe  our  recent 
success  in  finding  a  physical  foundation  for  thermal  fluctuations  that  yield  non-Nyquist  noise. 

2.  Non-Nyquist  Noise 

Nature  exhibits  several  types  of  noise  due  to  thermal  fluctuations  [18].  In  1827,  Brown  first 
reported  sporadic  motion  of  pollen  grains  in  water.  In  1905,  as  the  second  breakthrough  in  his 
Annus  Mirabilis,  Einstein  explained  this  “Brownian  motion”  by  assuming  random  impulses  from 
much  smaller  particles,  which  was  to  provide  the  first  definitive  evidence  for  atoms  and 
molecules.  As  a  function  of  frequency  (/)  Brownian  motion  exhibits  noise  with  a  power  spectral 
density  that  varies  as  S(f)oc  1  If  .  In  1926  Johnson  first  measured  electronic  noise  that  showed  a 
flat  spectral  density,  S(f)  =  const.  Nyquist  explained  this  “white”  noise  by  assuming  classical 
thermal  fluctuations  in  the  motion  of  the  electrons.  Also  in  1926  Johnson  measured  electronic 
noise  with  apparent  1  //  frequency  dependence.  Although  there  is  still  no  widely  accepted 
explanation  for  this  “ l/f  noise,”  empirically  it  is  the  most  common  low-frequency  behavior. 
Indeed  1// noise  is  found  in  metals,  semimetals,  semiconductors,  dielectrics,  ferroelectric s,  ionic 
conductors,  spin  glasses,  supercooled  liquids,  and  quantum  devices  [19-22],  as  well  as  in  music, 
speech,  neural  response,  and  human  perception  [23-26].  Although  no  single  model  is  likely  to 
explain  all  these  phenomena,  the  laws  of  nanothermodynamics  provide  a  unified  picture  for  1  // 
noise  in  many  materials.  Specifically,  the  general  principle  is  based  on  the  assumption  that 
particles  inside  local  regions  of  a  bulk  sample:  conserve  total  energy  by  including  non-extensive 
terms  in  E  (1st  law),  maintain  maximum  entropy  during  equilibrium  fluctuations  (2nd  law),  and/or 
exhibit  statistical  indistinguishability  of  identical  particles  consistent  with  quantum  mechanics, 
as  described  below. 

3.  Computer  Simulations  of  the  Ising  model 

Most  simulations  presented  here  utilize  the  Ising  model  for  binary  degrees  of  freedom 
(“spins”)  on  a  simple-cubic  lattice.  The  lattice  contains  a  total  of  N  spins,  with  exchange 
interaction  (7)  between  nearest-neighbor  spins,  and  periodic  boundary  conditions  between  the 
outside  surfaces  of  the  lattice.  Often  the  lattice  is  subdivided  into  smaller  regions  containing  n  < 
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N  spins  to  study  the  thermal  properties  of  small  systems  inside  a  bulk  sample.  Figure  2  (a)  shows 
the  net  magnetization  as  a  function  of  time  from  simulations  of  this  model.  Note  the  abrupt 
change  in  dynamics  at  time  t  =  0:  for  t  <  0  the  spin- flip  probability  is  governed  by  Boltzmann’s 
factor  alone  using  the  Metropolis  algorithm  Eq.  (1),  whereas  for  t  >  0  there  is  also  a  nonlinear 
correction  from  nanothermodynamics  Eq.  (2). 


e  AU,kliT  >  mnd[0,l) 

Boltzmann  Factor 

(1) 

,g(.s,„-s0)/kB  >  rcmd[0,  \) 

Nonlinear  Correction 

(2) 

In  Eq.  (1),  A U  is  the  change  in  interaction  energy,  ks  is  Boltzmann’s  constant,  and  T  is 
temperature.  In  Eq.  (2),  g  controls  the  strength  of  the  nonlinear  correction,  which  comes  from  the 
alignment  entropy  using  the  binomial  coefficient  for  binary  degrees  of  freedom:  Sm  /ks  = 
n\  . ,  „  .  n\ 


In 


[ '  ( n  +  m)\  ![4  ( n  —  m )] ! 


,  with  So  /k\i  =  In 


[(4«)!]2 


its  maximum  value.  When  g  =  0  (t< 0  in  Fig.  2 


(a))  systems  of  all  size  show  standard  Gaussian  fluctuations  characteristic  of  white  noise  (see 
Fig.  3).  When  g=l  (t> 0)  the  uppermost  set  of  data  (from  a  large  lattice  with  small  regions)  shows 
large  wandering  on  all  time  scales,  indicative  of  1// noise  (see  Fig.  3);  while  the  lower  two  sets  of 
data  (from  smaller  lattices  that  contain  a  single  region)  exhibit  sharp  jumps  and  steps 
characteristic  of  non-Gaussian  fluctuations.  Furthermore,  Fig.  2  (b)  shows  that  histograms  of  the 
simulations  exhibit  trimodal  behavior  (symbols),  similar  to  the  trimodal  behavior  found  from 
measurements  of  1// noise  in  a  spin  glass  and  ionic  conduction  through  a  nanopore,  shown  in  the 
upper  part  of  Fig.  2  (b)  (solid  lines).  In  contrast,  the  bottom  pair  of  lines,  from  fluctuations  in  two 
different  double- well  potentials,  shows  simpler  bimodal  behavior. 
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Figure  2.  From  Ref.  [17].  (a) 
Scaled  magnetization  versus 
time  at  two  temperatures,  from 
simulations  of  the  Ising  model. 
Note  that  simulations  from 
lattices  with  different  sizes,  N, 
are  offset  for  clarity.  The 
dynamics  changes  abruptly  at  t 
=  0  when  a  nonlinear  correction 
(Eq.  (2))  is  included  with 
Boltzmann’s  factor  (Eq.  (1)).  (b) 
Histograms  of  the  simulations 
from  N  =  12  (symbols)  and  from 
the  noise  measured  in  three 
types  of  samples  (solid  lines). 
Note  the  trimodal  behavior  from 
the  simulations,  spin  glass,  and 
nanopore  system. 
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Figure  3.  From  Ref.  [17].  Left  side  shows  frequency  dependence  of  noise  from  simulations  at  kBT/J=50  and 
500,  similar  to  those  in  Fig.  2.  S(f)  is  multiplied  by  N  to  scale  different  lattice  sizes  (given  in  the  legend)  and 


y. 


TIT. 


log (/)  is  multiplied  by  10  to  match  the  dB 
scale.  The  spectra  exhibiting  white  noise 
(bottom)  come  from  using  Eq.  (1)  alone. 
Spectra  exhibiting  l//-like  behavior 
(diagonal)  come  from  the  same  model 
using  both  Eq.  (1)  and  Eq.  (2).  Over  a 
broad  range  of  frequencies  these 
simulations  can  be  characterized  by  Sif)  oc 
Vf(T),  with  a(T)  ~  1.0  for  kBT/J=500 
(solid  line)  and  a(7)  ~  1.15  for  kBT/J =50 
(dotted  line).  Diamond-shaped  symbols, 
which  show  1  If  noise  at  low  frequencies 
and  white  noise  at  higher  frequencies, 
come  from  a  heterogeneous  system.  Right 
side  of  figure  shows  a(7)  from 
measurements  on  various  metallic  thin- 
films  (solid  symbols)  and  simulations 
(open  hexagons  connected  by  solid  lines). 
Note  that  the  temperatures  are  normalized 
by  T\,  where  a(7)  extrapolates  to  1. 


The  left  side  of  Fig.  3  shows  the  Fourier  transform  from  simulations  similar  to  those  in  Fig.  2 
(a).  Simulations  using  Eq.  (1)  alone  yield  white  noise  that  does  not  depend  on/(lower  symbols), 
whereas  adding  Eq.  (2)  yields  l//-like  noise  (symbols  along  the  diagonal).  In  fact  these 
simulations  can  be  characterized  by  S {f)oz  1  /faa\  with  a  temperature-dependent  spectral-density 
exponent  a(7)  that  is  consistent  with  the  measured  behavior  from  several  metals  at  lower 
temperatures,  as  shown  in  the  upper  corner  of  Fig.  3.  Similar  large-amplitude  low-frequency 
noise  is  found  in  most  substances.  Thus,  Eq.  (2)  provides  a  formula  for  obtaining  measured  1  //- 
like  noise,  including  deviations  from  pure  1 //"-behavior.  Moreover,  the  nonlinear  correction  is 
based  on  strict  adherence  to  the  fundamental  laws  of  thermodynamics,  as  described  below. 

Figure  4  depicts  some  of  the  mechanisms  that  justify  the  nonlinear  correction.  Consider  two 
binary  degrees  of  freedom,  e.g.  two  uniaxial  spins  that  can  be  up  or  down.  Figure  4  (a)  shows 
that  there  is  one  way  to  have  both  spins  up,  one  way  to  have  both  spins  down,  and  two  ways  to 
have  no  net  alignment.  The  alignment  entropy  is  obtained  from  Boltzmann’s  definition,  S/k b  = 
ln(f2),  using  the  multiplicity  of  each  alignment,  fL  The  dashed  (blue)  line  in  Fig.  4  (b)  shows 
how,  during  normal  thermal  fluctuations,  this  entropy  may  fluctuate  up-and-down  between  its 
maximum  S/k\>,  =  ln(2)  and  minimum  S/ks  =  0  values.  Although  seeming  to  violate  the  2nd  law  of 
thermodynamics,  various  explanations  have  been  proposed.  First:  the  entropy  of  small  system 
may  decrease  temporarily,  which  is  the  usual  assumption  of  various  fluctuation  theorems  [27]; 
but  clear  violations  of  the  2nd  law  require  completely-closed  systems,  which  cannot  be  measured. 
Second:  entropy  could  be  defined  using  Gibbs  ensembles  that  average  over  all  available  states, 
but  this  inhibits  using  entropy  for  time-dependent  and  out-of-equilibrium  behavior.  A  third 
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Figure  4.  From  Ref.  [17].  (a)  Sketch  of 
possible  states  in  a  region  containing 
two  binary  spins.  (a)  For 
distinguishable  spins  there  is  one  way 
to  have  both  up  (0+2=l)  or  down  (Q_ 
2=1),  but  two  ways  to  have  zero  net 
alignment  (O0=2).  (b)  During  thermal 
fluctuations  the  Boltzmann  entropy  of 
the  spins  goes  up  and  down  (dashed 
line).  To  maintain  maximum  entropy 
the  entropy  of  a  bath  must  go  down  and 
up  (dotted  line),  so  that  the  combined 
entropy  of  the  system  plus  bath  is 
constant  (solid  line),  (c)  When  the  bath 
has  high  entropy  each  low-entropy 
state  in  the  region  persists  twice  as 
long  as  expected  from  the  Boltzmann 
factor  alone,  (d)  Alternatively,  zero 
alignment  may  come  from  a  single 
state  that  contains  a  superposition  of 
spins,  consistent  with  de-localized 
particles  that  are  indistinguishable  in 
the  region. 


possibility  is  that  when  the  entropy  of  a  local  region  decreases,  the  entropy  of  its  bath  increases, 
so  that  the  total  entropy  remains  maximized.  Indeed,  Fig.  4  (b)  shows  how  the  entropy  of  the 
bath  (dotted  line)  may  balance  the  entropy  of  the  system  (dashed  line),  so  that  the  combined 
entropy  of  the  system  plus  bath  is  constant  (solid  line). 

This  entropy-transfer  mechanism  for  the  nonlinear  correction  comes  from  a  type  of  entropic 
force,  similar  to  Boltzmann’s  factor  [16].  Boltzmann’s  factor  favors  low-energy  states  because 
increasing  the  energy  of  the  system  lowers  the  entropy  of  the  bath.  Similarly,  the  nonlinear 
correction  favors  low-entropy  states  because  increasing  the  alignment  entropy  of  the  system 
lowers  the  entropy  of  the  bath.  In  fact,  for  the  two-spin  system  the  nonlinear  correction  precisely 
doubles  the  average  lifetime  of  each  fully  aligned  state,  as  shown  schematically  in  Fig.  4  (c). 
Thus,  each  aligned  state  becomes  as  likely  as  both  unaligned  states.  Information  theory  provides 
additional  insight  into  this  mechanism.  If  there  is  no  transfer  of  information  between  the  system 
and  its  environment,  then  the  alignment  of  the  system  cannot  be  known  and  its  multiplicity 
always  includes  all  four  states.  Again  the  entropy  is  constant,  remaining  at  the  value  Slk\>=\n{4). 

A  second  mechanism  for  the  nonlinear  correction  comes  from  the  statistics  of 
indistinguishable  particles,  as  shown  schematically  in  Fig.  4  (d).  To  match  the  probability  of 
each  alignment  in  Fig.  4  (c),  instead  of  doubling  the  likelihood  of  the  aligned  states,  the 
unaligned  state  could  be  a  single  superposition  of  up-  and  down-spins,  as  expected  for  identical 
particles  that  are  subject  to  the  exchange  interaction.  Indeed,  the  three  net  alignments  (up,  down, 
and  unaligned)  form  the  triplet  state  of  a  two-spin  system.  (The  less-likely  singlet  state  is 
neglected  in  such  a  simplistic  picture.)  Further  justification  for  this  interpretation  comes  from  the 
energy  and  its  fluctuations  shown  in  Fig.  5  (below),  where  a  similar  nonlinear  correction 
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minimizes  the  energy  of  small  regions  and  makes  their  entropy  extensive  and  additive,  consistent 
with  the  statistics  of  indistinguishable  particles.  Thus,  this  mechanism  for  extensive  entropy  in 
small  regions  is  similar  to  the  semi-classical  ideal  gas  that  resolves  Gibbs’  paradox  for  extensive 
entropy  in  large  volumes.  In  other  words,  the  nonlinear  correction  may  be  a  simplistic  way  to 
simulate  quantum-like  behavior  in  classical  models. 

A  third  way  to  understand  the  nonlinear  correction  is  from  conservation  of  energy  including 
Hill’s  subdivision  potential.  Consider  a  system  of  n  non-interacting  Ising-like  spins  with 
magnetic  moment  jub  in  field  B.  Each  spin  can  be  aligned  or  anti-aligned  with  B,  giving  energy 
-jUsB  or  +/i\>,B,  respectively.  The  single-particle  partition  function  is  Z\  =  e‘,H'llkBl  +e  _ 
Because  the  spins  are  non-interacting,  the  partition  function  for  n  spins  is  Z  =  (Zj)",  yielding  the 
free  energy  A  =  -nk^T  ln(Zi).  For  simplicity  let  B~$ 0,  so  that  A  =  -nksT  ln(2).  Again  using  the 
binomial  coefficient  for  the  exact  entropy  of  the  system,  the  internal  energy  becomes  Um=  A  + 

n\ 


TSm  =  ~nk\{T  ln(2)  +  k\{T  In 


[i(n  +  —  m)] 


.  If  the  system  is  at  its  average  alignment  m= 0, 


Stirling’s  approximation  for  the  factorials  yields  So  /kB  ~  n  ln(2)  and  Uo  ~  0.  However,  if  total 
energy  is  to  be  conserved  during  fluctuations,  there  is  a  non-extensive  contribution  to  internal 
energy  Um  =  Uo  -  Em,  where  Em  =  k\{T  (So  -  Sm)  is  Hill’s  subdivision  potential.  In  other  words, 
although  Uo  ~  0  at  m  =  0,  during  thermal  fluctuations  the  change  in  alignment  entropy  causes  a 
change  in  energy,  even  if  there  are  no  interactions  between  particles.  Fluctuations  in  Em  occur 
because  free  energy  is  a  thermal  average  over  all  states,  while  instantaneous  energy  and  entropy 
may  be  defined  for  every  state.  Note  that  when  m  ±  0,  Em  >  0  lowers  the  total  energy,  which 
favors  subdividing  a  bulk  sample  into  an  ensemble  of  regions  that  fluctuate  independently  to 
increase  the  fluctuations,  as  is  found  in  the  primary  response  of  most  materials.  Inserting  this  Em 
from  entropy  into  a  Boltzmann- like  factor  yields  the  nonlinear  correction,  Eq.  (2).  Interaction 
energies  that  appear  explicitly  in  Eq.  (1)  neglect  this  source  of  energy.  Thus,  Eq.  (2)  must  be 
included  if  total  energy  is  to  be  strictly  conserved  during  equilibrium  fluctuations. 

Deviations  from  thermodynamic  behavior  shown  by  standard  simulations  using  Eq.  (1)  alone 
are  due  to  finite-size  effects  from  assuming  homogeneous  systems  with  uniform  correlations. 
Specifically,  when  energy  fluctuations  are  averaged  over  a  volume  that  excludes  interacting 
particles  outside  the  volume,  correlations  across  the  interface  are  neglected.  For  large 
homogenous  systems  the  nonlinear  correction  gives  a  relatively  small  contribution  to 
conservation  of  energy:  n~>cc  yields  m-> 0,  and  Eq  ~  0.  However,  several  experimental 
techniques  have  established  that  the  primary  response  of  most  materials  comes  from  a 
heterogeneous  distribution  of  independently-relaxing  nanometer- sized  regions.  Indeed, 
dynamical  heterogeneity  in  the  correlations  between  regions  has  been  found  in  the  slow  response 
of  liquids,  glasses,  polymers,  and  crystals  [4-9].  In  fact,  extensive  studies  of  time-dependent 
specific  heat  at  low  temperatures  find  that  the  only  substance  to  show  no  evidence  for  localized 
excitations  is  chemically  and  isotopically  pure  NaF  crystals  [28].  Furthermore,  the  technique  of 
nonresonant  spectral  hole  burning  (NSHB)  establishes  that  excess  energy  is  localized  in  selected 
degrees  of  freedom  inside  a  sample  for  minutes,  or  even  hours  [29,30],  without  influencing  the 
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energy  in  neighboring  regions.  Thus,  complexity  in  the  primary  response  of  most  substances 
comes  from  thermodynamic  heterogeneity  due  to  an  ensemble  of  independently  relaxing  regions, 
consistent  with  nanothermodynamics,  but  different  from  the  behavior  assumed  for  standard 
thermodynamics. 

Energy  fluctuations  in  most  Monte  Carlo  simulations  exhibit  a  size  dependence  [31]  that 
differs  from  the  expectation  that  entropy  is  extensive  and  additive.  Figure  5  shows  the  average 
potential  energy  density  ( <U/J>/n )  and  its  fluctuations  (<(AU/kTy>/n)  as  a  function  of  the 
number  of  particles  ( n )  in  local  regions  of  a  large  lattice.  Again  the  simulations  utilize  the 
standard  Ising  model  for  binary  spins  on  a  simple-cubic  lattice  with  ferromagnetic  interaction 
between  nearest-neighbors.  The  solid  (black)  lines,  from  simulations  using  Boltzmann’s  factor 
alone  (Eq.  (1)),  show  constant  energy  density,  whereas  the  fluctuations  in  energy  density 
increase  with  increasing  region  size.  The  size  dependence  of  these  energy  fluctuations  yields  a 
size-dependent  specific  heat  that  is  inconsistent  with  NSHB  measurements,  while  integrating  the 
heat  capacity  yields  non-extensive  entropy  that  is  incompatible  with  standard  definitions.  The 
symbols  in  Fig.  5  show  the  energy  density  and  its  fluctuations  for  the  same  Ising  model  with 
various  strengths  ( g )  in  an  approximation  of  Eq.  (2)  that  yields  a  quadratic  correction  to 
Boltzmann’s  factor,  with  a  bypass. 

Quadratic  Correction  (3) 

The  bypass  comes  from  the  Kronecker  delta  that  gives  ( 1  -  SAU  0 )  =  0  when  A U  =  0.  A  practical 

reason  for  this  bypass  is  to  accelerate  slow  relaxation  and  avoid  frozen  response.  A  physical 
reason  is  that  if  A U  =  0,  spin  flips  can  occur  in  a  microcanonical  ensemble  without  coupling  to 
the  heat  bath.  Sufficient  disorder  in  a  local  region  may  remove  the  bypass  if  neighboring  states 
with  AU  =  0  are  not  available,  consistent  with  the  fact  that  1  If  noise  generally  increases  with 
increasing  disorder.  Furthermore,  most  materials  exhibit  a  combination  of  1  If  noise  at  low 
frequencies  and  white  noise  at  high  frequencies,  indicating  dynamics  without,  and  with,  the 
bypass,  respectively.  Indeed,  the  diamond- shaped  symbols  in  Fig.  3  come  from  a  heterogeneous 
mixture  of  Eq.  (2)  and  es'"^Sa)a^Smfi)'kB  >  rand[ 0,1).  In  any  case,  Figs.  5  (a)-(d)  show  that 
increasing  g  in  Eq.  (3)  reduces  the  energy  density  and  increases  the  energy  fluctuations  per 
particle,  until  g=l  where  the  energy  of  small  regions  is  minimized  and  the  specific  heat  is 
independent  of  region  size.  Thus,  g=l  yields  behavior  that  is  most  consistent  with 
thermodynamic  equilibrium:  energy  that  is  minimized  and  fluctuations  in  energy  that  yield 
extensive  entropy.  Empirical  evidence  for  extensive  entropy,  even  in  small  regions  inside  bulk 
samples,  comes  from  the  technique  of  NSHB  showing  that  the  specific  heat  is  constant  for 
independently  relaxing  regions  [32],  even  for  very  small  regions  having  n  ~  10  molecules  [33]. 
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Figure  5.  From  Ref.  [14].  Size  dependence  of  (a)  average  energy  per  particle  and  (b)  its  fluctuations  from 
MC  simulations  of  the  Ising  model  Solid  lines  come  from  using  the  standard  Metropolis  algorithm;  symbols 
come  from  simulations  with  varying  strength  of  a  nonlinear  correction,  with  g=  1  the  expected  correction  from 
a  Taylor  series  expansion  of  entropy.  The  dashed  lines  show  that  g=  1  yields  energy  reductions  proportional  to 
1  In  in  (a),  and  constant  energy  fluctuations  in  (b).  (c)  Average  energy  per  particle,  and  (d)  its  fluctuations,  as 
a  function  of  the  strength  of  the  nonlinear  correction.  Each  type  of  symbol  comes  from  a  different  region  size, 
as  given  in  the  legend. 
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4.  Comparison  of  Computer  Simulations  with  Measurements 

The  nonlinear  correction  also  improves  agreement  between  computer  simulations  and  the 
measured  response  of  many  materials.  The  symbols  in  Fig.  6  (a)  show  the  magnetic  susceptibility 
from  single-crystal  gadolinium  as  a  function  of  reduced  temperature.  At  high  temperatures  the 
data  have  a  slope  of  -1.0,  consistent  with  the  Curie-Weiss  law.  As  T~$TC  the  red  (dashed)  line 
approaches  a  constant  slope  of -1.24,  showing  that  standard  simulations  of  the  Ising  model  using 
Boltzmann’s  factor  alone  are  consistent  with  renormalization  group  theory  for  homogeneous 
samples  in  the  canonical  ensemble.  However,  close  inspection  of  the  data  reveals  that  the  slope  is 
not  constant  as  T~>Tc-  The  solid  (black)  curve  in  Fig.  6  (a),  which  shows  improved  agreement 
with  the  data  at  all  temperatures,  comes  from  simulations  of  the  same  Ising  model  with  the 
quadratic  correction  (g=l  in  Eq.  (3))  expected  for  heterogeneous  systems  obeying 
nanothermodynamics.  Figure  6  (b)  shows  the  difference  between  the  data  (symbols)  and 
simulations  of  the  Ising  model  without  (dashed  line),  and  with  (solid  line  at  the  origin),  the 
nonlinear  correction.  Indeed,  using  g=l  with  the  optimal  region  size  (n= 27)  reduces  the  standard 
deviation  by  an  order  of  magnitude.  Figure  6  (c)  shows  the  effective  scaling  exponent  yeff  as  a 
function  of  reduced  temperature,  from  the  magnitude  of  the  slope  when  plotted  as  in  Fig.  6  (a). 
The  data  (symbols)  and  simulations  with  nonlinear  correction  (solid  line)  show  conspicuous 
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Figure  6.  From  Ref.  [15].  (a)  Critical  scaling  plot  of  magnetic  susceptibility  versus  reduced  temperature,  r 
=  (T-Tc)/Tc.  Symbols  are  from  measurements  of  gadolinium.  Lines  are  from  simulations  of  the  Ising  model 
with  (g=  I ,  solid)  and  without  (g= 0,  dashed)  the  quadratic  correction  in  Eq.  (3).  (b)  Residual  plot  of  the  data 
and  g= 0  simulations,  from  g=  1  at  the  origin,  (c)  Effective  scaling  exponent  (magnitude  of  slope  from  Fig.  6 
(a))  for  the  Gd  data  and  Ising-model  simulations. 
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features  that  clearly  differ  from  the  monotonic  behavior  [34]  of  simulations  using  Boltzmann’s 
factor  alone  (dashed  line).  In  fact,  when  measured  to  within  0.01  %  of  Tc,  ultra-high-purity 
gadolinium  crystals  show  a  sharp  maximum  with  yeff  >  1.5,  and  yeff  1  as  T~>Tc  [35],  clearly 
incompatible  with  the  Ising  model  for  homogeneous  samples  in  the  canonical  ensemble. 

The  solid  line  in  Fig.  7  shows  the  excess  specific  heat  measured  in  the  colossal  magneto¬ 
resistance  material  LaMnCE.  The  peak  identifies  the  phase  transition  from  a  Jahn-Teller 
distortion  that  occurs  at  Tc  ~  730  K.  The  dashed  line  comes  from  simulations  of  the  Ising  model 
using  Boltzmann’s  factor  alone  (g=0),  while  the  solid  symbols  come  from  the  same  model  with 
the  nonlinear  correction  (g=l  in  Eq.  (3))  and  optimal  region  size  n=27.  Again  the  nonlinear 


correction  gives  significantly  better  agreement  with  the  measured  behavior. 


Figure  7.  From  Ref.  [16].  Excess 
specific  heat  as  a  function  of 
normalized  temperature.  Solid  line  is 
from  measurements  of  LaMnCL  near 
its  Jahn-Teller  distortion  temperature, 
Tc.  Dashed  line  and  symbols  come 
from  simulations  of  the  Ising  model 
without  (g=0)  and  with  (g=l)  the 
nonlinear  correction. 
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Figure  8.  From  Refs.  [16,37].  Pair  distribution  functions  (PDFs)  from  neutron  scattering  in  LaMnC^. 
Upper  lines  show  two  PDFs;  one  from  above  the  Jahn-Teller  distortion  temperature  and  one  from 
below,  with  an  identical  background  removed  from  both.  Lower  line  shows  the  difference  between 
these  two  PDFs.  Note  the  enhanced  correlation  at  short  distances,  then  abrupt  loss  of  correlation  at 
radius  r  ~  1.0  nm,  corresponding  to  the  distance  between  three  lattice  sites. 


Direct  evidence  for  heterogeneous  correlations  in  crystals  of  LaMn03  comes  from  neutron 
scattering  [36,37].  The  upper  two  lines  in  Fig.  8  show  the  pair  distribution  function  (PDF) 
measured  above,  and  below,  the  Jahn-Teller  transition  temperature.  The  difference  between  these 
PDFs  (lower  curve)  shows  strong  correlations  in  the  positions  of  neighboring  atoms  out  to  a 
distance  of  -1.0  nm,  then  abrupt  loss  in  correlation  beyond  [16]. 

Figure  9  shows  the  pair-correlation  function  from  simulations  of  the  Ising  model  without 
(g=0,  dashed  lines)  and  with  (g=l,  symbols)  the  nonlinear  correction.  When  y=0  there  is  a 
smooth  loss  in  correlation,  characteristic  of  the  classical  picture.  When  g=l  there  is  stronger 
correlation  over  the  three  contiguous  lattice  sites  within  each  region,  then  an  abrupt  loss  in 
correlation  to  the  neighboring  region,  consistent  with  the  measured  pair  distribution  function 
shown  in  the  lower  part  of  Fig.  8.  The  inset  of  Fig.  9  shows  that  linear  regression  on  three 
adjacent  sites  in  a  region,  from  simulations  with  g=l,  yields  a  correlation  length  that  is  similar  to 
the  radius  of  distinct  regions  deduced  from  neutron  scattering. 


r  (unit  cells) 


Figure  9.  From  Refs.  [16,37].  Pair 
correlation  function  versus  distance,  from 
simulations  of  the  Ising  model.  Dashed  lines 
show  smooth  loss  in  correlation  when  g= 0. 
When  g=1,  the  symbols  show  an  abrupt  loss 
in  correlation  between  regions,  with  a  more- 
gradual  loss  beyond  the  abrupt  jump.  Strong 
correlations  that  persist  for  three  lattice  sites 
are  consistent  with  Fig.  8,  and  n  =  33  in  Fig. 
7.  This  non-monotonic  loss  in  correlation 
deviates  from  the  classical  Ornstein- 
Zernicke  picture  [38],  where  the  pair- 
correlation  function  is  predicted  to  diminish 
smoothly  and  homogeneously  around  every 
atom.  Thus,  the  classical  picture  of 
homogeneous  correlations  is  invalid  on 
length  scales  of  nanometers,  where  quantum 
mechanics  may  influence  the  correlations. 
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Additional  direct  evidence  for  heterogeneous  correlations  comes  from  multi-dimensional 
NMR  [39,40].  Measurements  and  analysis  yield  an  average  size  for  the  independently  relaxing 
regions  of  10,  76,  and  390  molecules  (or  monomer  units)  in  glycerol,  ortho -terphenyl,  and 
polyvinyl  acetate,  respectively.  The  distribution  of  relaxation  times  gives  response  rates  that  can 
vary  by  several  orders  of  magnitude  across  nanometer  length  scales.  Nonresonant  spectral  hole 
burning  measurements  indicate  that  this  heterogeneity  corresponds  to  effective  local 
temperatures  that  also  vary  between  these  regions,  indicating  t/zermodynamic  heterogeneity. 
Thus,  accurate  models  and  theories  of  nanometer-sized  independently-relaxing  regions  must 
obey  the  laws  of  nanothermodynamics,  consistent  with  improved  agreement  between  computer 
simulations  and  measured  data  shown  in  Figs.  2,  3,  6,  and  7. 

5.  Simplified  Creutz  Model  for  a  Finite  System  with  Localized  Heat  Bath 

The  Creutz  model  consists  of  a  binary  system  coupled  to  an  explicit  heat  bath  [41].  The  model 
has  two  types  of  particles:  Ising-like  “spins”  that  have  fixed  positions  on  a  lattice  (the  system), 
and  kinetic-energy-carrying  “demons”  that  move  through  the  lattice  (the  heat  bath).  The  demons 
are  Einstein-like  oscillators  that  act  as  a  source  of  kinetic  energy,  with  energy-level  spacing 
equal  to  the  difference  between  Ising  states,  allowing  microcanonical  simulations  where  the  total 
energy  is  exactly  conserved  at  every  step.  The  dynamics  of  the  Creutz  model  involves  choosing  a 
demon,  moving  it  to  a  neighboring  site  in  a  chosen  direction,  and  testing  the  potential  energy  of 
the  spin  at  that  site  for  possible  inversion.  If  the  potential  energy  decreases  or  stays  the  same  the 
spin  always  inverts,  giving  any  extra  energy  to  the  demon.  If  the  potential  energy  tries  to 
increase,  the  demon  must  have  enough  energy  to  invert  the  spin  without  reaching  negative 
kinetic  energy.  Thus,  the  standard  Creutz  model  is  a  type  of  cellular  automata,  with  dynamics 
that  can  be  made  deterministic  [41].  We  have  investigated  a  simplified  Creutz  model,  where  the 
local  interaction  between  spins  is  replaced  by  a  uniform  interaction  with  an  external  field  (5) 
[42].  Key  advantages  of  the  simplified  Creutz  model  are  that  the  exact  canonical  distributions 
and  entropies  are  known  for  both  the  system  and  its  bath. 

The  zeroth  law  of  thermodynamics  states  that  when  two  systems  are  in  thermal  equilibrium 
they  have  the  same  temperature,  T,  but  the  law  does  not  specify  the  statistical  distribution  needed 
to  calculate  T.  Although  many  quantized  systems  can  be  characterized  by  Bose-Einstein  (BE)  or 
Fermi-Dirac  (FD)  statistics,  these  distributions  rely  on  the  usual  assumptions  of  Boltzmann’s 
factor  [43,15]:  the  heat  bath  must  be  effectively  infinite,  with  weak  coupling  between  the  system 
and  its  bath,  so  that  adding  one  unit  of  energy  to  the  system  does  not  change  the  probability  of 
adding  a  second  unit.  We  have  explored  some  consequences  of  using  a  finite  heat  bath,  with 
quantized  energies  for  the  system  and  bath.  We  find  that  assuming  the  BE  and  FD  distributions 
yields  unequal  canonical  temperatures  for  the  system  and  its  bath.  We  obtain  the  spin 
temperature  Ts  from  the  average  units  of  energy  per  spin  m/M,  and  the  demon  temperature  To 
from  the  average  units  of  energy  per  demon  n/N.  Specifically,  To  is  obtained  from  inverting  the 
BE  distribution,  e/ksTo  =  ln[(/z+A0/«],  while  Ts  is  obtained  from  inverting  the  FD  distribution, 
e/ksTs  =  In \(M-m)/m\.  From  computer  simulations  of  the  simplified  Creutz  model  we  find  that 


14 


To  and  T$  often  differ  significantly.  One  source  of  this  mismatch  comes  from  finite-size  effects, 
as  found  in  simulations  and  theory  of  critical  clusters  during  nucleation  [44,45].  In  the  limit  of 
large  systems  (N~M»  1 000)  this  temperature  mismatch  goes  to  zero  (To  ~^Ts  ),  as  the  energy 
densities  approach  the  canonical  BE  and  FD  distributions  n~>no  and  m-^mo.  However,  if 
energy-carrying  demons  are  reflected  internally  on  an  intermediate  length  (IL),  the  temperature 
mismatch  is  independent  of  system  size. 

We  model  this  persistent  mismatch  in  canonical  temperature  using  an  effective  interaction 
between  neighboring  quanta  of  energy,  combined  with  a  chemical  potential  to  accommodate 
conservation  of  energy.  Thus  this  temperature  mismatch  can  be  attributed  to  correlations 
between  units  of  energy  when  energy  is  quantized  and  conserved.  We  have  found  that  this  model 
adjusts  the  BE  and  FD  distributions  to  yield  a  single  T,  and  it  provides  good  agreement  with 
measured  excess  specific  heat  from  imperfect  crystals  at  low  T.  The  symbols  in  Fig.  10  come 
from  IF  dynamics  where  energetic  demons  reflect  at  boundaries  between  regions  with  sides  of 
length  L= 3  or  4,  as  given  in  the  figure.  Random  walk  dynamics  (RW,  not  shown  in  Fig.  10) 
would  yield  zero  offset  as  M~$  oo,  corresponding  to  the  dot-dashed  line  at  the  origin.  We 
emphasize  that  the  maximum  deviation  from  canonical  statistics  occurs  for  IF  dynamics  with 
L= 3,  returning  to  the  canonical  distributions  for  short -ranged  RW  or  for  long-ranged  ballistic 
motion.  This  IF  scale,  larger  than  atomic  spacing  but  smaller  than  sample  size,  may  be  related  to 
the  “universal”  mean- free  path  of  acoustic  waves  in  disordered  materials  [46].  Open  triangles 
show  the  offset  in  demon  energy  ( n-n0)/N  as  a  function  of  kBTD/s,  whereas  closed  triangles  show 
the  offset  in  spin  energy  (m-mo)/M  as  a  function  of  kBTs  /e.  Note  the  conspicuous  increase  in  (n- 
no)/N  and  shift  to  higher  kBTo  /£,  while  ( m-mo)/M  decreases  an  equal  amount  (to  conserve 
energy)  and  shifts  to  lower  kBTs/e.  Dotted  lines  connect  symbols  from  the  same  simulation  at  the 
same  total  energy.  Thus,  the  horizontal  offset  between  each  pair  of  connected  symbols  shows  the 
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Figure.  10.  From  reference  [42],  Difference 
between  simulated  and  theoretical  energy 
densities  as  a  function  of  canonical  temperature 
(triangles)  or  corrected  temperature  (circles). 
Open  and  solid  symbols  are  from  demons  and 
spins,  respectively,  with  symbol  shape  and  color 
identifying  the  distance  (L)  between  scattering 
surfaces  for  demon  motion.  Dotted  lines  connect 
pairs  of  points  from  the  same  simulation  at  the 
same  total  energy,  so  that  the  slope  of  these  lines 
shows  the  mismatch  in  canonical  temperature. 
Note  that  the  increase  in  demon  energy  above  the 
center  line  matches  the  decrease  in  spin  energy 
below  the  line,  as  required  by  conservation  of 
energy.  Other  lines  come  from  various  models,  as 
described  in  the  text.  Circles  show  the  corrected 
energy  differences  at  the  corrected  temperatures 
for  the  demons  (open  circles)  and  spins  (solid 
circles  near  the  center  of  the  open  circles).  The 
dot-dashed  line  shows  zero  deviation  from  the 
expected  energies. 
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mismatch  in  temperature,  which  exceeds  40%  for  L= 3  at  the  highest  energy.  Recent  results  [47] 
have  shown  that  in  fact  this  mismatch  in  temperature  comes  from  an  inhomogeneous  demon 
density  due  to  the  reflection  process.  Physically  the  inhomogeneous  particle  density  may  come 
from  standing  waves,  related  to  the  measured  “universal”  mean-free  path  of  acoustic  phonons.  In 
any  case,  the  maximum  mismatch  in  canonical  temperature  occurs  when  L= 3,  where  energetic 
demons  are  maximally  localized  to  the  single  site  at  the  center  of  each  3x3x3  region. 

The  solid  and  dashed  lines  in  Fig.  10  come  from  the  BE  and  FD  distributions  using  the  usual 
Creutz  temperature,  e/kBTo  =  ln[(no+A0//io].  Specifically,  the  upper  two  solid  lines  show  the  BE 
distribution  as  a  function  of  To,  1/  (ellk',T'>  -1)-  nJN,  which  agrees  with  simulations  of  ( n-rio)/N 
exactly  by  definition  of  To-  The  other  two  solid  lines  show  that  the  FD  distribution  as  a  function 
of  To,  1/  (V'//'"7''  +1)-  mJM,  has  a  small  increase  in  energy,  opposite  to  the  decrease  required  by 
conservation  of  energy  and  shown  by  the  simulations.  The  dashed  lines  that  roughly  follow  the 
lower  two  sets  of  symbols  again  use  To  in  the  FD  function,  but  with  the  numerator  and  energy 
spacing  adjusted  to  best  fit  the  data,  then  plotted  as  a  function  of  7s.  For  L= 3,  the  fit  parameter 
£=1.93+0.02  differs  from  the  actual  energy  spacing  a=  2,  while  the  numerator  0.835+0.005  also 
differs  significantly  from  its  actual  value  of  1.  Similarly,  some  values  of  this  FD  function  differ 
from  the  simulation  values  of  m/M  by  more  than  10  standard  deviations.  Thus,  for  IF  dynamics, 
the  energies  cannot  be  described  by  the  BE  and  FD  distributions  alone. 

To  characterize  deviations  from  the  canonical  distributions  we  use  an  interacting  lattice-gas 
model  for  adsorption  of  particles  on  neighboring  sites  [48],  adapted  to  treat  quanta  of  energy.  We 
find  that  the  offset  in  canonical  temperatures  can  be  modeled  by  an  effective  potential  w<0, 
which  causes  a  correlation  between  energy  quanta  on  neighboring  sites.  Although  agreement 
with  simulations  may  be  better  if  all  possible  arrangements  of  occupied  sites  in  a  region  are 
included,  here  for  simplicity  we  consider  only  nearest-neighbor  pairs.  The  partition  function  is:  f 
=  1  +  2 ae^tlkBTRW  +a2e^2elkBTRWe^w,l<BTRW  ,  where  a-e^^7™  is  the  energy-quantum  activity  at  the 
canonical  (random-walk)  temperature  Trw,  with  ft  the  chemical  potential.  The  three  terms  in  £ 
come  from:  no  units  of  energy,  one  unit  of  energy  on  either  site,  and  one  unit  of  energy  on  both 
sites,  which  includes  the  effective  interaction  potential.  The  probability  of  finding  no  energy  on 
either  site  is  1/f.  Using  conservation  of  total  energy,  the  density  of  energy  in  demons  and  spins 
becomes  n/N  =  mJN  +  Ale,  and  m/M  =  mJM  -  Ale,  respectively.  Here  we  use  A  as  an  adjustable 
parameter  that  governs  the  strength  of  the  energy  transfer  mechanism.  We  also  adjust  a  to  the 
best  value  for  each  set  of  simulations.  A  consistent  equilibrium  temperature  is  obtained  by 
inverting  either  energy  density,  with  the  appropriate  offset:  s/kBT  =  ln[  V(n/N-A/£  )+l]  or  e/kBT  = 
ln[l /(m/A+A/f )  -1].  The  open  circles  in  Fig.  10,  with  solid  circles  near  their  centers,  show  that 
the  corrected  BE  and  FD  distributions  give  good  agreement  with  the  energies  of  the  demons  and 
spins,  and  equilibrium  temperatures  that  match. 

Heat  capacity  in  the  standard  Debye  model  is  proportional  to  the  size  of  the  lattice,  so  that  at 
low  temperature  it  can  be  written  as  Cd  =  DMT.  Here  D  is  a  constant  that  depends  on  the 
number  of  degrees  of  freedom  per  lattice  site.  To  characterize  measured  excess  specific  heat  we 
assume  that  any  extra  energy  adsorbed  on  a  lattice  site  simply  adds  to  its  degrees  of  freedom,  so 
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that  D~>D+C(e /  £)(m  -  m0)/(MA).  Here  C  is  an  excess  specific  heat  constant,  (. m  -  mo)  gives  the 
additional  units  of  energy  at  each  site,  s  =  kgT 2  d\x\(cjldT  is  the  average  energy  from  the 

adsorption  model,  and  the  ratio  els  may  be  a  simplistic  way  to  count  the  extra  units  of  energy 
when  e  comes  from  a  distribution.  The  modified  and  normalized  specific  heat  can  be  written  as: 

Ca/(MT3)  =  D-C[2ae~e,kBT  +(2  +  »/£)fl2e'wd"'/v]/ £2  (4) 

Note  that  because  (2  +  w/s)  <  0,  the  quantity  in  square  brackets  is  often  negative,  yielding  excess 
specific  heat  for  C  >  0. 

Symbols  in  Fig.  11  show  the  measured  excess  specific  heat  reported  in  the  literature  for  three 
different  crystalline  substances.  Solid  curves  are  the  best  fits  to  these  data  using  Eq.  (4).  Despite 
the  simplistic  way  that  Einstein-like  oscillators  are  added  to  Debye-like  vibrations,  Eq.  (4)  gives 
good  agreement  with  the  excess  specific  heat  in  these  samples.  In  the  Creutz  model,  equal  energy 
spacing  simplifies  the  simulations  and  facilitates  conservation  of  energy.  From  the  measured 
excess  specific  heat  of  imperfect  crystals,  agreement  with  the  model  using  equal  energy  spacing 

connection  between  phonon  localization  and 

Figure.  11.  From  reference  [42], 
Normalized  specific  heat  as  a  function  of 
temperature  on  a  log-log  plot.  Symbols 
come  from  the  measured  specific  heat  of 
three  different  crystalline  substances: 
quartz  (black  squares),  ethanol  (red 
circles),  and  glycerol  (blue  triangles). 

Debye  theory  predicts  constant  values  of 
Cp/T3  at  these  low  temperatures,  as  seen 
below  3  K.  Solid  lines  come  from  best 
fits  to  the  excess  specific  heat  that  peaks 
above  10  K  using  Eq.  (1)  as  described  in 
the  text.  Agreement  with  the  data 
suggests  that  the  excess  specific  heat 
comes  from  extra  units  of  energy  that  are 
adsorbed  on  the  lattice  due  to  an  effective 
interaction  between  energy  units  on 
neighboring  sites. 


could  indicate  a  type  of  resonance,  or  other  inheren 
the  two-level  systems  [49] . 
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6.  Magnetoresistance  in  Nickel-Silicide  Nanowires  from  Nonequilibrium  Heating 

When  subjected  to  a  changing  magnetic  field,  the  low-temperature  resistance  of  Ni-Si 
nanowires  can  change  by  a  factor  of  two  or  more  [50-53].  Furthermore,  this  magnetoresistance 
shows  strong  hysteresis  as  a  function  of  field,  and  slow  relaxation  as  a  function  of  time.  We 
model  this  behavior  assuming  adiabatic  demagnetization  of  magnetic  spins,  with  charge  carriers 
that  are  also  weakly  coupled  to  the  thermal  bath. 

To  quantify  the  weak  thermal  link  we  calculate  the  local  temperature  and  its  change  relative  to 
the  bulk  temperature:  AT  =  TL  -  T .  For  adiabatic  spins,  thermodynamic  principles  [54]  predict 
that  changes  in  temperature  are  proportional  to  the  temperature-dependent  slope  of  the  magnetic 
moment  (dM/dT^n  divided  by  the  heat  capacity  at  constant  pressure  cp.  We  obtain  expressions 
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for  the  change  in  temperature  as  a  function  of  time  while  the  field  is  changing  t  <  tfm,  and  after 
the  field  sweep  is  stopped  t  >  tfin : 

AT  -  orrr(l-e~,/r)  t<tfin,  (5a) 

AT  =  arr{etfln'T  -l)^,/r  t  >  tfin.  (5b) 

Here  r  is  the  field  sweep  rate,  a  is  an  amplitude  factor,  and  r  is  the  time  constant  for  coupling  to 
the  thermal  bath.  Some  results  using  this  model  are  shown  in  Fig.  12. 
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Figure  12.  From  reference  [53].  Two-stage  magneto¬ 
resistance  relaxation  measurements.  Panels  (a),  (c)  & 
(e)  were  obtained  by  sweeping  the  field  from  zero  to 
+2  T,  and  then  holding  at  +2  T.  Subsequently,  panels 
(b),  (d)  &  (f)  were  obtained  after  allowing  the  sample 
to  equilibrate  for  a  long  time  (as  indicated  by  the 
offset  in  the  time  axis),  before  sweeping  the  field 
down  from  +2  T  to  zero.  In  each  of  the  panels  the 
gray  region  denotes  the  interval  over  which  the 
magnetic  field  was  swept.  Field  sweep  rates  are  0.2 
T/min  in  each  of  the  panels,  with  the  sample 
temperature  (T )  as  indicated.  Solid  curves  through 
the  data  represent  fits  from  Eqs.  5a  and  5b. 


7.  Deviations  from  the  Equipartition  Theorem  in  Simulations  of  Simple  Fluids 

Figure  13  shows  size-dependent  fluctuations  in  kinetic  energy,  normalized  by  a  formula  that 
would  yield  zero  if  the  system  obeyed  the  equipartition  theorem.  The  data  (symbols)  are  from 
microcanonical  molecular  dynamics  (MD)  simulations  of  systems  with  four  different  particle 
densities,  using  a  Lennard- Jones  (L-J)  interaction,  with  the  same  potential-energy  minimum  e  for 
all  systems.  The  symbols  come  from  averaging  the  behavior  at  three  temperatures  kT/e  =  0.01, 
0.02,  and  0.05,  with  error  bars  from  the  standard  deviation.  When  properly  normalized  by  kT/ne, 
these  fluctuations  exceed  the  equipartition  theorem  by  at  least  a  factor  of  ten.  For  all  systems,  the 
deviations  per  particle  decrease  with  increasing  system  size,  consistent  with  a  surface  term  that 
varies  as  n  In  (dashed  lines)  that  goes  to  zero  in  the  thermodynamic  limit.  However,  for  finite 
systems  these  size  dependent  fluctuations  cause  the  local  entropy  to  be  non-extensive  and  non¬ 
additive,  deviating  from  standard  thermodynamic  behavior,  and  disagreeing  with  measurements 
that  show  constant  specific  heat  for  independent  fluctuations  inside  bulk  materials. 
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Figure  13.  (unpublished)  Normalized 
fluctuations  in  kinetic  energy  as  a  function  of  the 
number  of  particles  in  the  region,  from  MD 
simulations  of  the  Lennard-Jones  model. 
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Deviations  from  standard  thermodynamics  shown  by  computer  simulations  arise  from 
assuming  homogeneous  correlations  between  particles.  Specifically,  when  energy  fluctuations 
are  averaged  over  a  volume  that  excludes  interacting  particles  outside  the  volume,  correlations 
that  influence  the  energy  fluctuations  are  neglected.  However,  several  experimental  techniques 
have  shown  that  primary  response  usually  comes  from  an  ensemble  of  independently-relaxing 
nanometer-sized  regions,  establishing  that  correlations  between  neighboring  particles  in  most 
materials  are  heterogeneous.  Indeed,  many  measurements  have  found  that  heterogeneity 
dominates  the  response  of  liquids,  glasses,  polymers,  and  crystals  [55-58].  In  fact,  detailed 
studies  of  time-dependent  specific  heat  at  low  temperatures  find  that  only  one  substance  agrees 
fully  with  the  Debye  T3  law:  chemically  and  isotopically  pure  NaF  crystals  [59].  Furthermore, 
results  from  the  technique  of  nonresonant  spectral  hole  burning  (NSHB)  indicate  that  specific 
heat  is  constant  for  independently  relaxing  regions  inside  bulk  samples  [60],  even  for  regions 
having  «<10  molecules  [61].  Thus,  the  primary  response  in  most  materials  is  dominated  by 
heterogeneous  correlations,  different  from  the  behavior  shown  by  standard  MD  simulations. 


8.  Discussion 

Some  remarkably  universal  empirical  formulas  are  often  used  to  characterize  the  measured 
response  from  many  materials.  For  example,  many  materials  that  exhibit  non-Nyquist  noise  have 
been  characterized  by  a  frequency  dependence  1  //“  with  a  ~  1,  as  shown  in  Figs.  2  and  3.  Non- 
classical  critical  scaling  has  been  used  for  behavior  near  phase  transitions,  as  shown  in  Fig.  6. 
Two-level  systems  of  uncertain  origin  have  been  used  for  the  excess  specific  heat  and  time- 
dependent  behavior  at  low  temperatures,  as  shown  in  Figs.  11-12.  In  many  cases  these 
expressions  are  convenient  mathematical  formulas  for  cataloging  the  measured  response  of 
complex  systems,  but  data  of  sufficient  quality  over  broad  enough  range  invariably  show 
deviations  from  these  formulas.  Many  models  have  been  proposed  for  each  of  the  empirical 
formulas,  so  that  the  deviations  may  be  a  decisive  way  to  distinguish  between  models.  The 
nonlinear  correction  to  Boltzmann’s  factor  provides  a  common  foundation  for  all  of  these 
formulas,  including  many  of  the  measured  deviations.  Moreover,  the  correction  is  necessary  to 
describe  the  thermal  equilibrium  of  any  sample  that  contains  independently  fluctuating  regions  if 
the  fluctuations  themselves  are  to  govern  the  distribution  of  region  sizes. 
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Standard  models  based  on  homogeneous  thermodynamics  have  been  unable  to  explain  several 
features  in  the  dynamics  of  complex  systems.  The  deviations  may  be  quite  subtle.  Indeed,  it  can 
be  difficult  to  see  curvature  in  the  data  on  a  log-log  critical-scaling  plot,  as  shown  in  Fig.  6  (a). 
Nevertheless,  other  researchers  have  also  recognized  that  most  ferromagnetic  materials  deviate 
from  standard  critical- scaling  behavior.  In  1989  Collins  wrote  [62]:  “The  critical  exponents  of 
iron  and  nickel  are  very  similar  to  each  other,  while  those  for  cobalt  are  clearly  different.  There 
is  no  theoretical  understanding  of  these  results.”  Also  in  1989  Hohenemser  et  al.  wrote  [63]:  “At 
the  same  time  our  review  makes  clear  that  when  one  restricts  the  analysis  to  the  best 
experiments,  only  a  few  materials  correspond  unambiguously  to  these  models,  while  most  do 
not.”  In  fact,  by  plotting  the  residuals  (as  in  Fig.  6  (b))  the  nonlinear  correction  clearly  improves 
the  agreement  between  the  Ising  model  and  measurements.  Moreover,  the  monotonic  behavior  of 
standard  simulations  of  the  Ising  model  cannot  match  the  sharp  temperature-dependent  features 
in  the  effective  scaling  exponent,  as  shown  in  Fig.  6  (c).  In  any  case,  such  nonlinear  corrections 
are  necessary  if  independently-fluctuating  nano  meter- sized  regions  inside  bulk  samples  are  to 
strictly  conserve  total  energy  and  maintain  maximum  entropy  during  equilibrium  fluctuations. 

9.  Conclusions 

We  have  shown  that  nonlinear  corrections  in  thermodynamics  provide  fundamental  insight 
into  the  local  dynamics  inside  complex  materials.  Although  full  understanding  of  all  features  will 
require  sophisticated  studies,  even  the  most-simplistic  models  can  match  measured  behavior 
when  appropriate  nonlinear  corrections  are  used,  as  shown  in  Figs.  2,  3,  6,  7,  9,  and  11.  From  the 
discussion  around  Fig.  4  it  is  possible  to  emphasize  the  following  results  regarding  the  laws  of 
nanothermodynamics.  The  1st  law  requires  that  total  energy  is  strictly  conserved  at  all  times, 
including  Hill’s  subdivision  potential  from  the  configurational  entropy  when  localized  regions 
fluctuate.  The  2nd  law  requires  that  the  entropy  of  an  isolated  system  must  never  decrease,  so  that 
if  a  local  region  fluctuates  into  a  low-entropy  state,  the  entropy  of  its  thermal  bath  must  increase 
to  compensate.  Thus,  the  nonlinear  corrections  extend  standard  thermodynamics  to  treat  finite- 
size  fluctuations,  while  strictly  maintaining  the  fundamental  laws  of  physics.  Moreover,  these 
corrections  yield  the  statistics  of  indistinguishable  particles  within  each  region,  as  is  needed  to 
avoid  non-extensive  entropy,  resolve  Gibbs’  paradox,  and  mimic  quantum- mechanical  behavior 
for  fundamental  particles  in  nanometer- sized  regions,  as  discussed  in  connection  with  Figs.  4  and 
5.  Indeed,  the  cooperative  dynamics  within  each  region  and  uncorrelated  dynamics  between 
neighboring  regions,  combined  with  the  fact  that  standard  thermodynamics  accurately  describes 
macroscopic  coherent  states,  suggests  a  connection  to  quantum  mechanics  [64,65].  Specifically, 
each  region  may  be  associated  with  localized  excitations  that  have  phase  incoherence  with 
excitations  in  neighboring  regions,  so  that  the  nonlinear  corrections  provide  a  fundamental 
connection  between  quantum  mechanics  on  the  scale  of  nanometers  and  the  bulk  behavior  of 
complex  materials.  Similar-  multi-scale  modeling  and  quantum-based  calculations  are  crucial  for 
accurate  simulations  of  energetic  materials  [66-68].  This  combination  of  large-scale  quantum 
behavior  with  small-system  thermodynamics  is  the  main  focus  of  our  current  research. 
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